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We demonstrate existence of non-pairwise interaction forces between vortices in multicomponent and layered 
superconducting systems. That is, in contrast to most common models, the interactions in a group of such 
vortices is not a universal superposition of Coulomb or Yukawa forces. Next we consider the properties of 
vortex clusters in Semi-Meissner state of type- 1.5 two-component superconductors. We show that under certain 
condition non-pairwise forces can contribute to formation of very complex vortex states in type- 1.5 regimes. 



I. INTRODUCTION 



The crucial importance of topological excitations in the 
physics of superconductivity made quantum vortex solutions 
in the Ginzburg-Landau (GL) models perhaps the most stud- 
ied examples of topological solitons (defined as localized 
lumps of energy characterized by a topological invariant) 1 ' 2 . 
These well studied vortex solutions are frequently used as 
generic testing objects for High Energy Physics and Cosmo- 
logical models 1 . In that broader context, especially spectac- 
ular theoretical works attempted on numerous occasions to 
identify topological solitons in the Skyrme and Faddeev mod- 
els with particles 1 . There the particle- solitons are lumps of 
energy which enjoy a topological protection against radiating 
their energy. Superconducting vortices, in spite being known 
to form a variety of "aggregate" vortex states: vortex liq- 
uids, glasses etc still do not show nearly as complex ground 
states as e.g. the Skyrme or Faddeev models 1 because of 
the smaller diversity of known intervortex interaction poten- 
tials. The recent works on multicomponent superconductors 
aimed at realizing more complex bound multi-vortex states in 
the so-called "type- 1.5 regime". In that regime, due to ex- 
istence of two components, thermodynamically stable vortex 
solutions were found that exhibit strongly non-monotonic in- 
teraction potentials between two vortices, with short-range re- 
pulsive and long-rage attractive parts 3 and thus form vortex 
clusters in low magnetic fields. After the recent experimen- 
tal publications 4 this physics received increased attention (see 
e.g. 6 ). 

In this work we investigate the structure of multi-vortex 
bound states in two-component superconductors beyond the 
validity of the linearized theory and find non-pairwise multi- 
body interactions which are especially pronounced at short 
intervortex separations. We show that in certain cases the 
non-pairwise forces are strong enough to led to extremely 
rich multivortex states. We also propose how these vortex 
states can be experimentally realized in layered superconduct- 
ing structures. 

We consider a Ginzburg-Landau model of a two-component 
superconductor which appears in various physical systems 
ranging from multiband and layered superconductors to pro- 
jected states of metallic hydrogen and models of neutron stars 
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Here e is the electric charge, ipi^ = 1^1,2 1^ 1,2 represent 
the superconducting components coupled by the gauge field 
A and the Josephson coupling rj. For a microscopic deriva- 
tion of such model for two-band superconductors see 7 . It was 
demonstrated recently in the framework or a self-consistent 
microscopic theory that this "minimalistic" two-band model 
describes qualitatively well intervortex interaction in multi- 
band superconductor in a rather wide range of temperatures 8 . 
We choose the units where H = c = 1. 

The key feature of the model is that one-flux quantum vor- 
tices induced by magnetic field are "composite". That is, 
they have a core around which both phases wind by 2tt, i.e. 
A0i } 2 = § c VOi^di = 2tt (C being a closed path around 
the vortex core). One-quanta vortex can thus be viewed here 
as a bound state of two "fractional" vortices 9-11 . The model 
exhibits type- 1.5 superconductivity when the magnetic field 
penetration length scale A is smaller than one of the char- 
acteristic length scales of the density variation £i 5 2 and also 
the conditions for short range repulsion and thermodynami- 
cal stability are satisfied 3 . In the type- 1.5 regime two vor- 
tices with similar circulation have interaction that is attrac- 
tive at long range (driven by density-density interaction) but 
repulsive at short range (due to current-current and magnetic 
interaction) 3 . This results in a formation of a "semi-Meissner" 
state in low magnetic fields where vortex clusters coexist with 
two-component Meissner state. The magnetization curves of 
a type- 1.5 superconductor are schematically shown on Fig. 1. 

We consider below two types of systems: (i) two-band 
and layered superconducting systems with weak and strong 
Josephson coupling and (ii) the systems where there is no 
Josephson coupling. The latter situation occurs in the con- 
text of the theories of liquid metallic hydrogen, deuterium 
and their mixtures 10 ' 11 , and physics of neutron star interiors, 
where the two fields represent electronic, protonic and E _ 
hyperon Cooper pair condensates 12 (for literature overview 
see 5 ). Because one cannot convert, say, electrons to protons 
or deuterons, the intercomponent Josephson coupling in this 
case is forbidden on symmetry grounds but the condensates 
are coupled by the vector potential. Experimentally the type- 
1.5-like physics can be also realized in a system of interlaced 
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Figure 1. A schematic picture of magnetization curves of type-I, 
type-II and type- 1.5 superconductors. Type- 1.5 superconductor has 
a first order phase transition in low magnetic fields. In low mag- 
netic field it can support a phase separation into vortex droplets and 
Meissner domains (the semi-Meissner state). 
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Figure 2. A schematic picture of a collection of Josephson-coupled 
layers of type-I and type-II superconductors. Vortices induced by 
magnetic field are kept in stacks by interlayer Josephson and elec- 
tromagnetic coupling. If in the type-I layers the cores will extend 
beyond the average magnetic field penetration length, these extended 
cores should cause attractive interaction between these vortex lines. 
The system can be used to model the type- 1.5 behavior. 



tex loop proliferation is based on London approximation 13 . In 
this approximation the fluctuations of densities are neglected 
and the intervortex interaction is then pairwise. Non-linear 
effects which lead to nonpairwise contributions in the inter- 
vortex interaction are, in general, less discussed. For cer- 
tain vortex configurations in a single-component Ginzburg- 
Landau model, nonlinear effects were studied recently in 14 . 

Here we demonstrate the importance of complicated non- 
pairwise forces between superconducting vortices arising in 
multicomponent systems. It has particularly important con- 
sequences for vortex clusters in the type- 1.5 regime (though 
also relevant for type-II two band superconductors). 



II. THREE-BODY INTERVORTEX FORCES 

First, let us present a highly accurate numerical study of 
a three body vortex problem. The interaction between vor- 
tices was investigated as follows (for a detailed description 
see apendix A): First a vortex pair is fixed in the center of the 
system. A third vortex is then inserted, and the energy is min- 
imized with respect to all the degrees of freedom, except the 
positions of the centers of the vortex cores in the first compo- 
nent. The procedure is then repeated for different positions of 
the third vortex. The resulting interaction energy is shown on 
Fig. 3 which shows pronounced violation of the superposition 
principle for intervortex forces in this system. To minimize 
the effects of discretization the calculation was performed on 
a ultra-high resolution grid of up to ex 10 7 points, with lattice 
spacing h ~ £i/100 (where £i is the coherence length of the 
dominant component in the limit of no coupling to the second 
component). We used 4-16 h on a 8-core cluster node to relax 
each data point in the interaction potential. We choose a "min- 
imally invasive" procedure of pinning a vortex on a numeri- 
cal grid, which gives most accurate long- and medium- range 
forces but, on the other hand, does not work for the (irrelevant 
for this study) very short intervortex distances (for details see 
Appendix A). Consequently in Fig. 3 no data is given for too 
closely placed vortices. 

In all the type- 1.5 regimes which we considered we found 
diverse and pronounced to various degrees non-pairwise inter- 
actions. 



III. VORTEX CLUSTERS IN A SEMI-MEISSNER STATE 
AND NON-PAIRWISE INTERACTIONS. 



Josephson-coupled type-I/type-II multilayers shown schemat- 
ically in Fig. 2. 

The usual framework within which vortex matter in super- 
conductors and superfluids is usually discussed relies on the 
assumption that interactions in a system of vortices is a su- 
perposition of pairwise forces. Indeed the most usual analy- 
sis of interaction between well separated topological solitons 
involves linearization. By nature of this approximation, the 
interaction in a system of multiple vortices is a superposition 
of two-body forces. Similarly, for example the description of 
phase transitions in type-II superconductors in terms of vor- 



Let us now investigate how the presence of non-pairwise in- 
teractions along with non-monotonic two-body forces affects 
the magnetic response. Below we report highly accurate nu- 
merical solutions for N- vortex bound states in several regimes 
(for technical details see Appendix B). Animations showing 
the evolution of the system, during the numerical energy min- 
imization, from the various initial configurations to the vortex 
clusters are also available as supplementary material 15 . 

The figures (Fig. 4-Fig. 11) show various bound states of 
multiple flux quanta in U(l) x U(l) as well as in Josephson- 
coupled U{1) models. Consider first the case where < 0. 
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Figure 3. Interaction energy between a single vortex and a fixed vortex pair (position of the fixed vortices are shown by two black lines). GL 
parameters are : (ai,/5i) = (-1.0, 1.0) and e = 1.3; (a 2 ,/3 2 ) = (3.0, 0.5), rj = 3 for panels a, b , c ; (a 2 , AO = (3.0, 0.5), r\ = 7 for d , e , 
f panels and (a 2 , #2) = (—0.0625, 0.25), 77 = 0.5 for panels g , h , i . Top row (a , d , g ) gives the total interaction energy, while second one 
(b , e , h ) gives the sum of pairwise interactions. The third row (c , f , i ) displays the difference between these energies, which demonstrates 
the existence of strong nontrivial, non-pairwise, (three body) interactions. The regimes shown in the left, middle and right columns are type-II, 
type- 1.5 and type-I correspondingly. 



There appear very interesting geometrical properties of the 
vortex ground states (shown on Fig. 4,Fig. 5). One can clearly 
see that with growing number of vortices, the local vortex 
structure strongly depends on the number of vortices in a clus- 
ter. The striking feature which is to various degrees is manifest 
on all the figures is the coexistence and competition between 
type-II-like behavior of the first condensates which attempts 
to form a regular vortex lattice and type-I-like behavior of the 
second condensate. Namely the second condensate mimics 
the formation of a single large normal domain. Also like in 
a genuine type-I superconductor, this component has super- 
current density which is predominantly concentrated on the 
boundary of the domain. Clearly it energetically prefers a for- 
mation of a circular boundary. These competing tendencies 
in the type- 1.5 system result in neither hexagonal nor circular 
boundary. The next visually striking effect is that the vortex 
solutions are very different inside the vortex cluster and on 
the cluster boundary. This shows up especially clearly in the 
current density. 

The multibody forces in this nonlinear theory, originate in 
nontrivial deformation of vortices by their neighbors in a vor- 
tex cluster. First the qualitatively new physics which arises in 
these clusters is the appearance of gradients of the phase dif- 
ference V(#i — 62) between the condensate fields. It is clearly 
seen on the panels f from the plotted quantity Im^*^)- One 
of the mechanisms of the generation of the phase difference 
which we observed was associated with splitting of the vortex 



cores of the components ^2 driven by competing interac- 
tions. It leads to a "dipole"-like configuration of the phase dif- 
ference of the two components which in turn results in contri- 
butions to muti-body forces associated with the induced phase 
difference gradients. This splitting exists in cases of zero, as 
well as finite Josephson coupling, though it is smaller in the 
later case. We observed also more complicated configurations 
like "quadrupoles" of the phase difference fields. These con- 
figurations occurred when the vortices broke their axial sym- 
metry by relegating more phase gradients in the areas where 
density was depleted by neighboring vortices. Another source 
of multibody forces was associated with nontrivial condensate 
densities modulations for a group of several vortices. 

Note that the presence of gradients in the phase difference, 
along with gradients of the relative density of two conden- 
sates is known to lead to contributions from self-generated 
Skyrme-like terms to magnetic energy density 16 . This makes 
the physics of the two-component vortex cluster boundary and 
the resulting multibody forces, in general, a complicated non- 
linear problem. 

With increased number of vortices the simulations fre- 
quently produced long-living bound states of irregular shapes 
like that shown on Figs. 10 and 11. 

Parameter sets where the system is close to type-I regime 
(like the one shown of Fig. (7)) show that the transition (in the 
parameter space of the model) to type-I regime manifests itself 
in a depletion of current densities inside vortex clusters and 
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Figure 4. The ground state of a N v = 9 flux quanta configuration in type-1.5 £7(1) x £7(1) superconductor (i.e. r\ — 0). The parameters 
of the potential are (ai, = (—1.00, 1.00) and («2, /fe) = (—0.60, 1.00), while the electric charge is e = 1.48. The physical quantities 
displayed here are a the magnetic flux density, b (resp. c) is the density of the first (resp. second) condensate \ipi } 2 \ 2 • d (resp. e) shows the 
norm of the supercurrent in the first (resp. second) component. The panel f is Im^i ^2) which is nonzero when there appears a difference 
between the phases of ip* and ^2. The second component has a type-I like behavior: its density is depleted in the vortex cluster and its current 
is mostly concentrated on the boundary of the cluster. 
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Figure 5. The ground state of a N v = 12 flux quanta configuration. The parameter set is the same as in Fig. 4. Here the global 8-folded 
discrete symmetry of the cluster has been broken down to the 3-folded symmetry. There is a competition between type-I-like (normal circular 
cluster with a boundary current) and type-II-like tendencies (vortex lattice). 



relative increase of current density on the cluster's boundary. 



In type-I regime the system forms one circular cluster where 
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Figure 6. Ground state of 9 vortices in a U(l) x U(l) superconductor with two active bands {0.2^2) — (—0.6, 1) (without interband 
coupling). The parameter set here is as in Fig. 4 except e — 1.55 which places the system closer to the transition with the type-I regime. This 
is manifest in the more circular boundary of the cluster compared to Fig. 4. 
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Figure 7. Ground state of an N v = 9 vortex configuration with the parameter set given by Fig. 4, but with e = 1.59 and added Josephson 
coupling 77 = 0.1. Although the Josephson term introduced an energy penalty for phase difference, it has little effect on the vortex clus- 
ter boundary. At the boundary the competing magnetic and density-density interactions win over phase-locking terms and generates phase 
difference gradients. Also in this case the system is close to type-I regime: i.e. most of the current is concentrated on the boundary of the 
cluster. 



the current is concentrated on the boundary (i.e. a single giant vortex). 





Figure 8. Ground state of 9 vortices in a superconductor with two active bands. Parameters of the interacting potential are (ai,/?i) = 
(—1.00, 1.00), (c*2, P2) — (—0.0625, 0.25) while the interband coupling is r\ — 0.5 which is substantially larger than on Fig. 7. The electric 
charge, parameterizing the penetration depth of the magnetic field, is e = 1.30 so that the well in the nonmonotonic interacting potential is 
very small. In this case there is visible admixure of the current in second component in vortices inside the cluster. 
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Figure 9. Ground state of N v = 12 vortices, for a superconductor with broken U(l) x U(l) symmetry, with the parameters as in Fig. 8. The 
system compromises between both type-I-like and type-II-like tendencies in the different superconducting components. 



Note that in single-component type-I superconductors the domains of normal phase have a circular boundary only when 
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Figure 10. Elongated cluster of vortices in the same superconductor as in Fig. 8, but with N v = 18 vortices. 




Figure 11. Bound state of N v = 28 vortices, for a superconductor with broken U(l) x U(l) symmetry, with the parameters as in Fig. 8. This 
irregularly shaped cluster represents a local minimum of the free energy. The local minima originate from the competing interactions yielding 
a complicated free energy landscape. 



the effects of stray fields outside sample are neglected. In fields typically lead to formation of macroscopically large 
finite samples of type-I superconductors the stray magnetic stripes of the normal phase rather than circular domains 2 , 
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though other geometries were also observed 17 . Similarly in 
realistic experimental setups especially for type- 1.5 bilayers, 
stray magnetic fields could lead to vortex stripes rather than 
circular vortex clusters formation. 



IV. NON-COMPACT VORTEX CLUSTERS AND 
NON-PAIRWISE INTERACTIONS 

Next we study the structure formation in a regime with rela- 
tively strong non-pairwise forces. In this section we consider a 
situation where the passive (i.e. with positive o^) second band 
is coupled to the first band by extremely strong Josephson 
coupling r] = 7.0 (shown on Fig. 12). This coupling imposes 
a strong energy penalty both for disparities of the conden- 
sates variations and for phase difference. The two and three 
body forces in that regime are shown in the middle column 
on Fig. 3. Indeed the two-body forces in this case have only 
a weak non-monotonicity. Importantly, there is an anisotropy 
of three body forces shown in Fig. 3. It clearly diminishes the 
energetic benefits of a triangle-like states compared to line- 
like vortex states. Identifying the ground states in this regime 
numerically (which involves four-body and higher order inter- 
actions) is much more complicated than in the previous cases. 
We get a flat and complicated energy landscape and the out- 
come of the energy minimization strongly depends on initial 
configuration (see Appendix for the description of numerical 
procedures). The Fig. 12 shows the typical non-universal out- 
come of the energy minimization in this case. Striking feature 
here is formation of vortex stripe-like configuration. Indeed it 
strongly contradicts the ground state structure predicted by the 
two-body forces in this system. Namely the axially symmet- 
ric two-body potentials with long range attraction and short- 
range repulsion (which we have in this case) do not allow 
stripe formation in the ground state configurations.. There- 
fore the observation of the vortex stripes signals that the struc- 
ture formation, along with short-range character of attractive 
tail, is influenced by repulsive multi-body forces in these cases 
(in contrast to the structure formation in the previous section 
which was dominated by two-body forces). Note that even in 
this regime, the system exhibits self-induced gradients of the 
phase difference, in spite of the strong Josephson coupling. 
In order to study the role of initial conditions we consider the 
vortex ordering in a similar regime but starting with 30 vor- 
tices placed at distances larger than the minimum of the two- 
body potential. If there were only two-body forces this vortex 
clusters would contract to minimize the energy. In contrast in 
the energy minimization process the vortex cluster first ex- 
pands (see the animation of the evolution of the system in the 
numerical energy minimization process in the supplementary 
material 15 ). Subsequently it breaks into a few sub-clusters and 
vortex chains. At the final stage of the evolution each cluster 
contracts. Final intervortex distances in each sub-cluster is 
smaller than intervortex distance in the initial state shown on 
Fig. 13. 

Formation of highly disordered states and vortex chains due 
to short-range nature of the attractive potentials and many- 
body forces was a generic outcome of the simulation in the 



similar type- 1.5 regimes with strong Josephson coupling, in 
spite of negligible effects of ultra-fine numerical grid. 

In this section we considered the regimes where the attrac- 
tive part of two-body interaction was relatively weak and of 
short range. Physically, the fact that there are also multi-body 
forces which energetically penalize the hexagonal arrange- 
ment in large groups of vortices implies that at finite tempera- 
tures small clusters and irregular chain-like structures should 
easily form as well for entropic and kinetic reasons. 



V CONCLUSIONS. 

In conclusion, in this work we investigated structures of 
vortex clusters in the Semi-Meissner state and demonstrated 
that non-pairwise interaction forces can be especially impor- 
tant in multicomponent and layered superconductors. It re- 
sults in the very rich physics associated with the vortex clus- 
ters in type- 1.5 regime. Namely N-quanta clusters can be 
quite different from simple superpositions of N single vortex 
solutions. In general they should have multiple local minima 
in the energy landscape. Thus vortex clusters can have ex- 
tremely irregular shapes even in absence of vortex pinning. 
The discussed above regimes can be directly probed experi- 
mentally in type-I/type-II multilayers where intercomponent 
coupling can be tuned and be distinguished from pinning ef- 
fects by IV characteristics. The found here non-pairwise 
forces should similarly be important in type-II multiband su- 
perconductors or layered structures for understanding com- 
pact configurations of pinned vortex clusters. 

We also remark that in certain cases the separation into 
vortex and Meissner domains also implies phase separation 
into domains with different broken symmetries. Consider 
U(l) x U(l) model. Even if the second component there 
is not completely depleted in the vortex cluster, its density 
is suppressed and as a consequence the magnetic binding en- 
ergy between vortices with different phase windings (A#i = 
2tt, A<9 2 = 0) and (A6 1 = 0, Af9 2 = 2tt) can be arbitrarily 
small 9 . Moreover the vortex ordering energies in the compo- 
nent with more depleted density will also be small. As a re- 
sult, even small thermal fluctuation can drive vortex sublattice 
melting transition 10 in a macroscopically large vortex droplet. 
In that case the fractional vortices in weaker component tear 
themselves off the fractional vortices in strong component and 
form a disordered state. Note that the vortex sublattice melt- 
ing is associated with the phase transition from C/(l)xE/(l)to 
U(l) state 10 . I.e. that vortex cluster where one sublattice has 
melted will represent a domain of U(l) phase (associated with 
the superconducting state of strong component) immersed in 
domain of vortexless U(l) x U(l) Meissner state. 

We thank A. Gurevich and J.M. Speight for stimulating 
discussions. The work is supported by the NSF CAREER 
Award No. DMR-0955902, by Swedish Research Council, 
and by the Knut and Alice Wallenberg Foundation through 
the Royal Swedish Academy of Sciences fellowship. Compu- 
tations were performed at the National Supercomputer Center 
(NSC) in Linkoping, Sweden. 
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Figure 12. A bound state of an N v = 25 vortex configuration in case when superconductivity in the second band is due to interband proximity 
effect and the Josephson coupling is very strong 77 = 7.0. The initial configuration in this simulation was a giant vortex. Other parameters are 
(ai,ft) (-1.00, 1.00), (a 2 , A2) (3.00, 0.50), e 1.30. 
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Figure 13. A bound state of an N v = 30 vortex configuration for the system with parameters like in Fig. 12 which was obtained using a dilute 
initial configuration of vortices 15 . 



Appendix A: Finite difference energy minimization terns and inter- vortex interaction energies are found by mini- 

mizing this functional subject to relevant constraints, such as 
vortex positions. To do this numerically, we discretize the 
To calculate intervortex interaction energies we use finite system on a regular grid. To have the numerical results unbi- 
difference energy minimization. Ground states of vortex sys- ased we use a non-adaptive grid where the grid spacing h is 
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the same everywhere in the domain. The Hamiltonian is then 
discretized using the finite difference approach: 
Gradients are defined as 



(V/) M +i = 



/(* + !)-/(*) 
h 



and magnetic flux is computed by line integration 



B, 



h 2 



A • dr 



(Al) 



(A2) 



where u is the square with the corners z, i + 1, j, j + 1. The 
energy density in the grid point (i, j) then depends on function 
values in z, j and its neighbors. 

The optimization scheme which is used in the first part 
of the paper is a modified version of the Newton-Raphson 
method. To minimize boundary effects we use free boundary 
conditions. Vortices are inserted using various initial configu- 
rations. 

In order to calculate the inter- vortex interaction energy, we 
have to fix the position of vortices. Fixing a vortex position 
requires a special care to avoid the situation where the pinning 
substantially affects the vortex solution. We fix the vortex po- 
sition by the following method. In the vortex center the con- 
densate density is zero. We then fix the density only of the 
central point of the dominant component | xpi | of the vortex to 
be zero in a given position of numerical grid. This effectively 
prevents the vortex from moving but does not prevent core 
splitting of |^i | and \ip2 \ due to competing interactions. This 
"point pinning" method also has advantage of being a "min- 
imially invasive" since only the position of the core singular- 
ity is fixed. Thus it allows calculate medium- and long-range 
forces with greatest accuracy. 

However, at the same time, obviously, this method does 
not work for too short intervortex separation. For too short 
vortex separation it leads to the following easily identifiable 
artifact: a vortex core of one of the vortices elongates to be 
zero at both pinning centers allowing the second vortex to 
unpin and escape, while satisfying the energy minimization 
constraint. This behavior can be easily remedied by differ- 
ent pinning schemes. However for consistency and also be- 
cause very short-distance intervortex forces are irrelevant for 
the questions studied in this paper we use one pinning proce- 
dure. 

Convergence is determined as follows: 

1) We choose a particular grid spacing hi and number of 
grid points Ni = Ni x ■ Ni y giving a system size L x = h ■ 
(N\ x — 1), Ly = h ■ (Ni y — 1). Then we minimize the 
energy until it does not change in a few thousand iterations. 
This gives us E(h\). 

2) We decrease grid spacing h by a factor of 2 or 3 while 
retaining the system size L x , L y using spline interpolation. 
Then, we once again iterate until the energy does not change 
in a few thousand iterations, giving us E(h2) and so forth. We 
then determine convergence from 



We use grid sizes up to N « 10 T which gives very high accu- 
racy, typically C < 10 -4 . 



Appendix B: Finite element energy minimization 



In the second part of the paper we use the uncontrained 
energy minimization. Bound state vortex configurations are 
minima of Ginzburg-Landau energy (1). This means that 
functional minimization of (1), from an appropriate initial 
state describing several flux quanta, should lead to bound 
state (if it exists). We consider the two-dimensional problem 
T = j Q F defined on the bounded domain Q C R 2 , supple- 
mented by 'open' boundary conditions. 

The variational problem is defined for numerical pur- 
pose using a finite element formulation provided by the 
Freefem++ 18 framework. Discretization within finite element 
formulation is done via a (homogeneous) triangulation on f2, 
based on Delaunay-Voronoi algorithm. Functions are decom- 
posed on a continuous piecewise quadratic basis on each tri- 
angle. 

Contrary to the numerical method used in the first part, the 
accuracy does not depend only on the 'number of grid points' . 
The accuracy of such method is controlled through the num- 
ber of triangles, (we typically used 10 5 ), the order of expan- 
sion of the basis on each triangle (P2 elements are 2nd order 
polynomial basis on each triangle), and also the order of the 
quadrature formula for the integral on the triangles. 

Once the problem is mathematically well posed, a numeri- 
cal optimization algorithm is used to solve the variational non- 
linear problem (i.e. to find the minima of J 7 ). We used here 
Nonlinear Conjugate Gradient method Algorithm is iterated 
until relative variation of the norm of the gradient of the func- 
tional T with respect to all degrees of freedom is less than 
10 -6 . To be sure that our results are not numerical artifacts 
of this particular minimization scheme, we also performed 
standard Steepest Descent calculations and checked it leads 
to similar results. 

Minimization starts with an initial guess: a field configura- 
tion carrying the N v flux quanta described by 



*" = u *f[Jl ( 1 + tanh (^(^(*> y) - tS) ) e 



A — — (sin 6, — cos 6) , 



(Bl) 



E(h n ) - E(h n+ i) 
E(h n ) 



= C. 



(A3) 



where a = 1 , 2 , u a is the vacuum expectation value of each 
scalar field, the parameter £ a give the core sizes while 9 and 
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1Z are 



N v 



&i(x,y) = tan" 



i y-Vi 



N v 



Ki(x,y) = ^(x-Xi) 2 + {y-y i )' 2 



(B2) 



(xi,yi) are the initial position of a given vortex. Then, all 
degrees of freedom are relaxed simultaneously without any 
constraint to obtain high accurate solutions of the Ginzburg- 
Landau equations. 



The initial guess (Bl) allows starting from various very dif- 
ferent initial configurations, depending on the values of the 
fa, yi). Since we know from two-body calculations, the pre- 
ferred distance between two vortices, we choose to start ei- 
ther in the repulsive or in the attractive tail of the two-body 
interaction potential. Animations showing the evolution of 
the system from these various initial configurations, during 
the described above energy minimization is available as on- 
line supplementary material 15 . Note that we do not solve a 
dynamical problem here, the main purpose of these movies 
is that they reflect the structure of the free energy landscape 
and to give intuitive illustration of intervortex forces. For 
a reader interested in the a relationship between this energy 
minimization dynamics and the dynamics of Time Dependent 
Ginzburg-Landau model we refer to Ref. 19 . 

The method described in the first part of the paper was also 
used for unconstrained simulations to doublecheck some of 
the results. 
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